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Abstract 

'a': 

Using pole decompositions as starting points, the one parameter (— 1 < c < 1) nonlocal and 
nonlinear Zhdanov- Trubnikov (ZT) equation for the steady shapes of premixed gaseous flames is 
studied in the large- wrinkle limit. The singular integral equations for pole densities are closely related 
to those satisfied by the spectral density in the O(n) matrix model, with n = — 2 They can be 
solved via the introduction of complex resolvents and the use of complex analysis. We retrieve results 
obtained recently for — 1 < c < 0, and we explain and cure their pathologies when they are continued 
naively to < c < 1. Moreover, for any — 1 < c < 1, we derive closed-form expressions for the shapes 
of steady isolated flame crests, and then bicoalesced periodic fronts. These theoretical results fully 
agree with numerical resolutions. Open problems are evoked. 

PACS: 47.70.Pq ; 47.20.Ky ; 05.20.Jj ; 02.30.Em. 
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1 Introduction 

Most flames propagating into premixed gaseous reactants may be viewed as fronts: their actual 
thickness I ~ V -Dti^ch ~ Ui^ch (based on flat-flame speed ul, heat diffusivity Z?th, and chemical time 
t c h) is often much smaller than the wavelength and amplitude of their deformations. Being also very 
subsonic (i.e. ul <C speed of sound), such combustion fronts border fluids of constant and uniform 
densities: p u on the fresh side, and pb < p u in the burnt gas. As the Atwood number A = — jj^ < 1 is 
nonzero, such flames are subject to the hydrodynamic, hence nonlocal, Darrieus |Dar) -Landau |Lan44) 
(DL) wrinkling instability at large wavelengths. At shorter scales though, the variations in local 
normal burning speed u n (relative to reactants) with front mean curvature [Mar51| bring about a 
neutral wavelength L n eutrai proportional to £, yet much longer in practice IPC82] . Viewing wrinkled 
flames as fronts embedded in incompressible flows largely facilitates numerical simulations CM11 
and theoretical analyses thereof. 

Acknowledging that A -C 1 implies a weak DL instability - and hence asymptotically small front 
slopes and slow evolutions - Sivashinsky analysis |Siv77] provided the first systematic weakly-nonlinear 
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description of the local amplitude of a wrinkling: a companion numerical work |MS77j confirmed that 
the Michelson-Sivashinsky (MS) equation (Eq. ((TJ) with c = 0) correctly captures the slow spontaneous 
dynamics of flat-on-average flames if A — > + . In physical situations, A ranges from 0.2 — 0.5 for 
thermonuclear flames in Supernovae [RHN04 to 0.65 — 0.85 for chemical fronts. It is thus important 
to incorporate higher orders in A. The first two subleading orders can be absorbed in a renormalization 
of coefficients of the MS equation [SC8 7, KL02, Kaz05 . At the third one a new quadratic term enters 
the equation, but the fourth order correction can again be absorbed in coefficients [KPPj . A similar 
nonlinearity first appeared in derivations of an equation for the wrinkling amplitude that assumed 
A = 0(1) and postulated a small flame slope |ZT89[ IBKL99] . 

If a single space coordinate is retained, the Zhdanov- Trubnikov (ZT) type of equation [ZT89 so 
obtained for slow spontaneous evolutions of front wrinkles reads as follows, in rescaled form: 



where c is some function of A, 4>(t, x) is a dimensionless front deformation about a flat shape (0 = 0); 
the subscripts denote derivatives with respect to rescaled time (t) or coordinate (x). For space- periodic 
cells, < v < 1 will denote the neutral-to-actual wavelength ratio, provided the reference length is 
suitably chosen. For the isolated crests of infinite wavelengths, v > could be rescaled to unity by 
a change of variables but it is kept in ([T]) for comparisons with cells. The dependence of u n — ul on 
curvature gave the v$ xx term in ([1]), while the Hilbert transform: 



in the right-hand side encodes the nonlocal DL instability: a small <f>{t,x) oc exp(tnt + inx) indeed 
has w = \k\ — vk 2 as growth rate. The stabilizing local nonlinearity in (|T|) combines geometry [Siv77 
and hydrodynamics |SC87j . The nonlocal one, absent from the MS equation (c = 0), is mainly 
fluid-mechanical since 4> 2 + Ti.[4> x ] 2 is proportional to the kinetic energy of wrinkling-induced flow 
disturbances. It may be destabilizing if c < 0, or overstabilizing if c > 0, affecting the solutions of ([TJ) 
in several regards, e.g. at the front tips (local, usually sharp maxima of </>). 

Besides an interest per se, solving the nonlinear and nonlocal ([1} analytically with c as general as 
possible - and viewed here as a free parameter - might help one fit flame shapes from experiments 

QSDG012 ; or from simulations that use A = O(l) [CMllj . This could also form the basis of 
stability analyses of curved fronts, e.g. generalizing that of |VM00j to large wrinkles. The task is 
facilitated by the fact that, just like the MS equation [TFH85 , ([TJ where x is replaced by a complex 
variable Z = x + iB, admits meromorphic solutions with pairs of complex-conjugate simple poles 
Z_k{t) — Z£(t). When restricted to the real line Z — x, these solutions are real- valued and represent 
the front slope. Then, ([1]) is converted into coupled nonlinear differential equations for the poles Z^it) 
|Jou91j . When the front acquires a steady shape, those poles of <j) x that remain at rest for t — » oo 
generically are aligned along parallels to the imaginary axis, due to the same mechanism as for c = 

TFH85] : at Zk = i-Bfc (mod 2tt) when x-periodicity is assumed, and also at Z m = it + ib m (mod 2tt) 
if Neumann boundary conditions are employed (or periodic ones with x 4-> —x symmetry). In any 
such steady configurations, the equations that govern the 2N motionless pole altitudes may be viewed 
as expressing an electrostatic (or kinematic) equilibrium under mutual pairwise repulsions between 
charges (or fluid sources), all subject to an external field originating from the DL instability mechanism. 
When N 3> 1, corresponding to large wrinkles, the spacing between consecutive poles scales like 
Bn/N, enabling one to adopt a continuous approximation of the equilibrium conditions, and rewrite 
them as integral equations for pole densities [TFH85I . As first shown for c = in [JD08 , and confirmed 
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when c < [JD12] , a key step to describe the structure of large steady wrinkles is to first solve the 
simpler case of isolated crests ; these have infinite wavelengths and poles at = i-Bfe in arbitrary 
number 27V. 

In |JD12j . such integral equations were solved when — I < c < by exploiting peculiar properties 
of index- i Meixner-Pollaczek orthogonal polynomials evidenced in |Dun85j . firstly for isolated crests 
and then for monocoalesced periodic fronts (a single maximum of 4> per wavelength). In both cases, 
the results nicely agreed with numerical resolutions of the discrete pole equations when TV 3> 1. But 
the naive continuation to c > led to unsatisfactory, yet intriguing results. For example, the front 
slope so obtained for isolated crests with c > reads 
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4>x(x) = psin 7sinh ( — , (3) 



with 73* = ^ vN and c = tanh 2 (7T7/2) ; this lies surprisingly close to the numerically-determined 
slope profile whenever < c < 0.6, except very close to x = where ((3]) predicts an essential 
singularity. The pole density associated with (|3]) also stays very close to its numerical counterpart, 
except at \B\ <C 73*, where it oscillated infinitely many times and becomes negative at places: this is 
disallowed, for pole densities must be nonnegative. As conjectured in |JD12j . such localized pathologies 
might partly result from the failure of an assumption implicitly made in the continuation from c < 
to c > 0, namely that the density support would still include 73 = 0, instead of being limited to 
< 73 m ; n < |73| < 73 max w 73* for some 73 m j n > 0, and one may imagine that 73 m i n /73 max increases 
rapidly with c. An analysis of the discrete pole equations for c — > 1~ supports the first conjecture, 
and the results of Sec. |3]|4] will confirm and quantify both. 

To this end, another method is devised to compute the flame shapes analytically. The key step 
is to introduce an appropriate resolvent, which can be interpreted as a complex electric (or velocity) 
field caused by charges (or fluid sources) distributed according to the unknown pole density profile. 
The equations that rule the positions of the poles then turn into a nonlocal, scalar Riemann-Hilbert 
problem for the resolvent, which might be solved by complex analysis techniques. The advantage of 
this method is that, once the resolvent is found, the pole density and the front slope can be obtained 
easily, without extra integration. 

In this way, we reproduce the results obtained in |JD12j for — I < c < 0, and explain and fix the 
pathologies of (J3j) for c > 0. The method is valid whatever the sign of c, and generalizes results on 
wrinkled fronts with two unlike crests per wavelength, which was so far available for the MS equation 
only |JD08j . In all the configurations envisaged, we compare the analytical predictions with numerical 
determinations of the discrete T3& and of </>. 

The paper is organized as follows. In Sec.[2j the pole decompositions and their continuous approx- 
imations are introduced. We introduce at a formal level the notion of resolvent and we give equivalent 
forms of the ZT equation in § 12.31 The resolvent approach is then exploited for isolated crests (Sec. [3]), 
and next adapted to two-crested rc-periodic fronts (Sec. H} ; in both cases comparisons with numerical 
results are presented and discussed. We summarize our results and formulate open problems concern- 
ing flame shapes in Section [5] A resolution of the discrete pole equations for c — > 1~, and technical 
details about resolvent determinations, are presented in Appendices lAl and [El respectively. 

We also find links between flames fronts governed by the ZT equation, and other statistical- 
mechanical problems inspired from random matrix theory. Since the main body of the article is 
devoted to flames, we postpone this discussion to Appendix C. 
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2 Pole decompositions, isolated crests and the limit of large 
crests 



2.1 Discrete poles 

Once properly shifted as to be centered, steady isolated crests of infinite wavelengths have <f> — <p(x) 
where 4>x is a sum of 27V simple pole contributions oc x _\ Bk , \k\ = 1, . . . , N. B-f. = —Bk ensures that 
4>(x) is real valued when x is real. Since H[ x l iB ] = ^zrg > the dominant balance in ([T]) implies that 
each pole has residue jz^'- 

o N 1 

-2v \ ~\ 1 



— y — — . (4) 



fc=-jV K 
fe^O 

For ([I]) to be satisfied, the pole altitudes Bk obey coupled equations, originally expressed for k = 
±1,...,±N as |Jou91j : 

^ 1 - csgn(g fc )sgn(gj) = sgn(g fc ) 

.f^ (l-c^Bk-Bj) 2v ■ U 

Without loss of information, ([S]) may be rewritten in terms of Bx, . . . , B^ > only. For k = 1, . . . , N: 

N 1 t , JV 

El 1 + c 11 

whereby the discontinuous sgn(B) function no longer shows up. 

Summing B fe -(Eq. ©) yields < Bi H + B N = vN{^ - 1), whence c < 1. This must 

be considered jointly with the constraint — 1 < c needed from B\jv > (obtained for c — > — 1, cf. 
[JD12I ). For any fixed — 1 < c < 1, the -Bfc's will scale like vN if N — > oo, as in the c = case, 
whereby the typical spacing Bk — Bk-i is ~ v ~ B^/N. One may thus introduce a density 

e (B) > (7) 

here such that g(B)dB measures the number of poles with positive altitudes in [B,B + dB]. This 
continuous approximation enables to rewrite © as an integral equation, valid for B 6 [-B m in, -B max ]: 

B — gCgQdg; _ n f B — e (B')dB' = 1_ 

Bmin B — B' 2 J Bra . a B + B' 2v ' U 
where j denotes the Cauchy principal value, and the constant n is defined as: 

1 

n=-2 r 



(9) 



Notice that — 1 < c < is mapped to — 2 < n < 0, and it will be convenient to set in this regime 
c = — tan 2 (7T7/2), so that: 

0<7<l/2, n = -2cos(tt 7 ). (10) 

When < c < 1, we rather have n < —2, and a convenient parametrization is c = tanh 2 (717/2), so 
that: 

< 7, n = -2cosh(7T7). (11) 
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The equation <j8j) must be solved subjected to the normalization condition: 



/ 



B, 



g{B) dB = N, 



(12) 



and the continuous version of Q will express the front slope 4> x in terms of g(B) as: 



( 



/ 

J B, 



B, 



g{B)dB 



) 



<t>x(x) 



1 -c 



Im 



\x - B 



(13) 



The above g{B) may be nonzero for B > only, contrary to the two-sided even density p(B) = 
g{B) + g(—B) employed in |JD12| . Also notice that ([5]) allows for a lower end -B m ; n > of the density 
support, besides the upper end -B m ax ; both have to be determined as parts of the solution. For c < 0, 
B mm happened to vanish |JD12j . whereas an argument presented in Appendix [X] shows that B nlln > 
in the limit c — > 1~ . One of the outcome of this article is that B m - ln > for any < c < 1 (see Fig.Q}, 
while B m i a = for any — 1 < c < 0. 

2.2 Preliminary comments 

The kind of systems ([5|) , which express a repulsion between the poles Bk and some of their images 
(here the B_fc = — Bk) with possibly different intensity (here, if k, k' > 1, the interaction of Bk and 
B-k' is (— n/2) times stronger than the mutual interaction between Bk and Bk>), is ubiquitous in 
statistical physics. In the context of quantum integrable system, they typically appear as the result of 
a Bethe Ansatz to determine the eigenvalues of transfer matrices. This repulsion phenomenon is also 
observed for the zeros of orthogonal polynomials [MMFMG07 , optimal approximation points and the 
eigenvalues of random matrices Dei98, ForlO . The connection between those topics and flame fronts 
is not accidental: we point out in Appendix [C] that the ZT equation coincides with the large N limit 
of a Schwinger-Dyson equation in a statistic- mechanical ensemble of N particles called " /3-deformed 
O(n) matrix model". 

The problem of describing how poles satisfying an equation like (JS|) condense when N is large, i.e. 
of solving singular integral equations like ([8]) , is sometimes more convenient to consider from the point 
of view of complex analysis. There often exists an appropriate Riemann surface containing the support 
of g, on which the density can be continued to an analytic function g. One can then use the powerful 
tools of algebraic geometry to determine it, in particular the behavior of g at its singularities fixes 
much of the solution. We may bring up in the same spirit the techniques of conformal mapping used to 
solve Dirichlet problems for harmonic functions, which is well-known in electrostatics. The techniques 
to solve (J5|) effectively have been developed prior to this work in the context of matrix models, for all 
values of c [GK89 ( Kos89, EK95 ( EK96 . In this article, we adapt them to the determination of flame 
shapes. 



Before going into details, it is noteworthy that the behavior of g at the edges of the supports 
features some universality: it depends only on the type of short-distance repulsion between poles, and 
does not on the details of external potential (the right-hand side in ©). If Bk remains away from 
their images, the density will typically behave as a squareroot near B m i n > 0. With (|8]) this happens 
for c > 0, and obviously the flame slope 4> x (x) will be smooth near x = 0. But if B m i n — 0, the poles 
close to are subjected to a strong self-interaction but also a strong interaction with their images at 
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_2l I i I i I i I i I i I i I I 

-6 -4 -2 2 4 6 

x coordinate 

Figure 1: Comparison of theoretical and numerical pole density (top, theoretically given by (|43jl ) and 
flame slope (bottom, theoretically given by ([IB))) for an isolated crest with c > 0. The plot assumes 
c = 0.5, N = 100 and 1/u = 199.5. 
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opposite position since the latter are also close to 0. The exponent a such that g(B) oc B a when 
B — > + now depends continuously on the relative strength of interactions (see Sec [ 



I 2 

a = 7 = ~cos -1 (-n/2) = - tan -1 (V^) (14) 

7T 7T 

With (|SJ), this results for — 1 < c < in a cusp of the flame profile at x = 0, such that (f> x (x) oc 
— sgn(:r)|a;| -7 . 

Anticipating on Sec. we find that the exact solution of ((SJ) nicely agrees with what the numerical 
resolution of the discrete pole equations (j6|) gives for N ^> 1, and with the slope 4> x (x) ensuing from 
(|13jl . In the regime — 1 < c < 0, this exact solution is constructed in § 13.31 and was already known 
and compared to numerics in [JD121 Fig. 3] ; the corresponding flame profile has a cusp as we just 
discussed. For < c < 1, the exact solution is devised in § 13.21 and displayed in Fig. [TJ Unlike jJD121 
Fig. 5] where the pathological formula ((3|) was plotted, the theoretical prediction in Fig. [1] does not 
feature a singularity at x — 0, as it should be (we used in both plots the exact same parameters). 

2.3 Continuum of poles and resolvents 

A function <p x (x) (assumed square- integrable) can always be represented uniquely as |Tit37j : 

<t>x(x) = T^ImW(k) (15) 
1 — c 

where W(z) is a holomorphic function in the domain {Rez < 0}, called resolvent. In this paragraph - 
which can be read independently - we clarify what does it mean for <p x (x) to solve the ZT equation in 
terms of the resolvent, and therefore put the discrete pole decompositions Q in a broader perspective. 
A more elementary approach will be adopted in § 13.11 

Let us look for solutions of ZT such that W(z) is actually holomorphic in C\C, where C is a reunion 
of bounded arcs in the region {Rez > 0}, which is invariant by complex conjugation. We allow C 
to touch finitely many times the axis iR: these points of contact would correspond to singularities of 
<fi x (x). Besides, we assume that the front is asymptotically flat, i.e. W(z) = o(l) when z — > oo. Then, 
we can use Cauchy residue formula to represent: 

(16) 



c 



where we introduced the density: 

= W(B-i0)-W(B-i0) 
2m 

The discrete pole decomposition ((4]) corresponds to the special case where g(B) is a sum of N Dirac 
masses located at B = Bk for k = 1, . . . , N. Here, we rather want to consider the case where g(B) is 
continuous. 

The assumption C — C* allows to rewrite (fT5")) as: 

<fe (z) = — [W(ix) - W(-ix)] , (18) 

when x € M is not at a singularity of <p x {x). Thanks to the assumed analytic properties of W(z), the 
Hilbert transform of (j) x (x) can be computed explicitly: 

U[ct> x ]{x) = [Wiix) + W(-ix)] . (19) 
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Therefore, the ZT equation is equivalent to: for any z € iM, 

W 2 (z) + W 2 {-z) + nW{z)W{-z) + W'(z) + W'(-z) - (l/v)[W(z) + W(-z)] = 0. (20) 

where n was defined in j9]). And, since (|20[) is an equality between analytic functions, it must be valid 
in the whole domain of analyticity, namely for any z € C \ C. 

Large wrinkles correspond to the limit v — > but with W(z) scaling like 1/v so as to keep a 
macroscopic front slope in (fT5"|) . In this limit, the derivative term in (|2"U|) can be neglected and we find 
£ \yV\{z) — 0, where we set: 

£[W](z) = W 2 (z) + W 2 (-z) + nW(z)W{-z) - (1/v) [W(z) + W(-z)] . (21) 

Computing the discontinuity of £[W](z) at a point B G C, we find: 

[W(B + iO) - W{B - i0)] x [W(B + iO) + W{B - iO) + nW{-B) -l/v] = 0. 

Therefore: 

W(B + iO) + W(B-iO)+nW(-B) = 1/v. (22) 

In terms of the density g(B), this equation is equivalent to (JSJ which has been derived in ^ 12.11 for the 
continuum limit of a discrete pole configuration. 

Conversely, imagine that we have a function W(z) is holomorphic in a domain C \ C, which 
decays at infinity, which has a continuous density g(B), and which satisfies ([22j) . It follows that 
£[M /r ](z) is holomorphic on C \ (C U — C) and decays at infinity. Furthermore, by multiplying (|2"2"j) by 
W(B + iO) — W(B — iO) and rearranging, we find that £ \W\{z) is continuous across C. By parity, it 
is also continuous across — C, and therefore, £ [W](z) is an entire function decaying at infinity. By 
Liouville theorem, it must vanish: in other words, (f> x (x) given by (1151) is solution to the ZT equation 
([1]) in the limit of large wrinkles. 



3 Resolvents and isolated crests 
3.1 Strategy 

Even without refering to § 12.31 the form of © suggests to introduce the so-called resolvent, defined 
for a complex variable z as: 



W(z) = J^ ^± B ~. (23) 

In the context of an electrostatic (of a fluid-mechanical) analogy, W(z) would represent the complex 
electric (or velocity) field generated by charges (or fluid sources) deposited according to the density 
g(B) along the segment [B m - ln , B mayi } in the z-complex plane. Equivalently, W(z) is the unique function 
which is holomorphic on C \ [-B m ; n , B max ], behaves as N/z when z — > oo, and is discontinuous on 
[-^min , Bmax] with a given jump: 

WjB-jV-WjB + m = 
2m 

The Cauchy principal value in ((5J) can then be expressed as: 

g(B') dB' W{B - iO) + W(B + iO) 
/; . W = 2 ■ (25) 
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We can also reformulate ([8]) in terms of the resolvent only: for B £ [B m i n , B max ], 

W(B + iO) + W(B - iO) + nW(-B) = 1 /v. (26) 
The normalization condition (| 12|) is translated into the requirement: 

W(z)~N/z, z->oo, (27) 
and the front slope can be expressed readily by comparing of (fT5|) and (I23[) : 

^x(x) = ^Im[W(ia;)]. (28) 

The strategy to solve (f2l)| has been developed in |EK 95 , EK96] . The first step is to assume that 
the position of the density support [-B m in, -B ma x] is given, and determine it by consistency only at the 
end. The second step consists in performing a change of variable which uniformizes the complex plane 
with two cuts [B m i n ,B max ] and [— -B max , — -B min ], since points on these two segments are involved in 
(1261) . In other words, one constructs: 

= W(z(^)), (29) 

where ip belongs to some domain T> in the complex plane, so that (|26p relates boundary values of ui 
on T>, in such a way that its resolution becomes "easy" using Schwarz reflection principle. This can 
be done with trigonometric functions if -B m m = 0, whereas elliptic functions shows up if B m i n > 0. 
The third step is the resolution of the equation for w taking into account its analytical properties and 
the normalization (fT2j). which result in the determination the position of the support [B m i n , B max ] . 
If several solutions for the support are available, there is in general a unique one which ensures 
positivity of the density ((TJ. We already announce the result that B m i n = whenever — 1 < c < 0, 
while i? m i n > in the case < c < 1, and we now explain the implementation of this method in both 
cases. The limit cases c = —1, 0, 1 are better discussed separately, respectively in |JD12] . |JD08j and 
Appendix [2] 

3.2 Solution for < c < 1 

We assume in this paragraph -B m i n > 0, and later find this implies < c < 1. It is convenient to 
choose as change of variable: 

z = B min snk(ip), fc = min , (30) 

^max 

where sn^ is the Jacobi elliptic sine function [GW89] . Let K and K' be the complete elliptic integrals 
with modulus k. The z-complex plane with its two cuts [— B max , — -B mm ] U [-Bmin, -B m ax] is mapped onto 
the rectangle T> of vertices ±K and ±iK' (see Fig. [2]). More precisely: the segment [-B m i n > -Bmax] ± iO 
is mapped to the half-edge [K, K ± \K'] in the rectangle ; the lower (resp. upper) side of the segment 
] — oo, —-Bmax] U [B max , +oo[ is mapped to the upper (resp. lower) horizontal edge of the rectangle, 
and in particular the point z — oo correspond to tp = ±iK' ; and the segment ] — -B ma x, — -B mra [±i0 is 
mapped to [-K, -K ± iK'}. 



u(ip) = W(z(ip)) defines a holomorphic function in the rectangle T>. The analytical properties 
of the resolvent (|2l?l) provide relations between values of lu on the boundary of T>. Firstly, since the 
resolvent is continuous across ] — oo, — B max ] U [B msK + oo[, we must have: 

u){i>) = ui(ip + 2iK'), (31) 



9 



2-pLanc 



-K + iK'_ 
-if 



-if - iif' 



iK' 



-iff' 



6 
K 



Figure 2: Mapping between the z-plane (top panel) with its two cuts on [— B mull — B max ] U [B m - ln , 5 m ax] 
and the rectangle T> in the ^i-plane (lower panel). For special points in the -0-plane, we indicate below 
(in red) the corresponding value of z. The dashed blue and dotted green paths indicates how the 
mapping "opens" the cuts. 
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when ifj belongs to the lower horizontal edges of T>. Similarly, the continuity of W(z) across the 
segment [-B max , -B min ] implies: 

ui(ip) = u{-2K - i>), (32) 
when i/j belongs to the edge [-K — \K' 1 —K + iK']. Eventually, (j2"6")l itself is equivalent to: 

u]{f) + u){^-AK) + xiu(ii-2K) = l/v, (33) 

when ip belongs to the edge [K — iK', K + iK']. 

The key observation is that successive applications of Schwarz reflection principle and ([3T j) - ([33"l) 
allows to continue analytically u> on the whole -0-complex plane, so that (|3ip ~(|33 p are satisfied ev- 
erywhere. Indeed, (f3"Tj) allows to define u>(ip) as a 2ii4T'-periodic function on the strip S = {—K < 
Reip < K}, while (|32l allows to define uj(ip) on the doubled strip S 2 = {— 3K < Reip < K}. Then, 
(1331) provides a relation between values of lo on the boundary of S. Furthermore, if we take a point 
ip € 1S2 — the points ip — 2K and ip belong to 6>2- So, enforcing (|33l) defines as an analytic 

function of ip € 5^ = S2 U (S2 — which satisfies (|33)) everywhere on S' 2 . This argument can be 
used recursively to analytically continue u(tp) on all strips of the form £2 + AmK for any integer to, 
so as to cover the whole complex plane. 

The problem is now reduced to an easier one, namely finding the most general analytic functions 
which satisfy (|3"T j) - (j33]) . i.e. with an appropriate behavior under 2K and 2ii^' translations of their 
argument. The general solution can be put in the form: 

v{2 + n) 

where g(ip) is an arbitrary analytic function which is 2K and 2iK' periodic, and uj+(ip) is any analytic 
function satisfying: 

uj+{^ + 2iK') = w+(il>), lu+{4j + 2K) =e^LJ+{ip), (35) 

Here, the parametrization pT|) n = — (e^ 7 + e^ 1 ) arises quite naturally. In the constant term 

is an obvious particular solution of the inhomogeneous equation (|33[) . while the second term is the 
general solution of (|31l) - (|33|) with vanishing right-hand side. 

Notice that, since W was holomorphic on CU{oo}\ [B m i n , B max ], and the right-hand side of (|33p is 
a constant, uj(ip) must be an entire function. Therefore, in the solution we are looking for, g(^)Lj + (i/;) 
(which also satisfies (l3"5j) ) is entire as well. We explain in Appendix IB1 that this is impossible, unless: 

jK — 2pK for some positive integer p. (36) 

In particular, 7 must be real, hence c > 0. Besides, when (I36[) holds, we easily guess a holomorphic 
solution of (|55j) : 

u + (%l>) = e***l K '. (37) 

With this choice, we are looking for an entire function gfy) which is biperiodic, hence bounded, hence 
constant by Liouville theorem. We thus obtain 

w(V0 = J +n) + G cosh[7r( 7 /2 + i>/K% (38) 

for some constant G. The normalization (|27p implies two consistency equations given by: 

N 

= + - (ip - iK') + o(tp - iK '). (39) 

ip^nK' i>max 
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Jointly with (|36|) . they provide three equations for the three unknowns B m i D , S max and G. The integer 
p seems to allow for a discrete set of solutions, but we now argue that p = 1 is the only one which 
is physically admissible. Indeed, let us come back to the density of poles (|2"4"|) . After our change of 
variables, it can be expressed as: 

e m) = u(K-ix(B)K')-Lj(K + i x (B)K') 

2lTT 

where < x(B) < 1 is the unique point for which ip(B) = K + \x{B)K' satisfies ([50)1 for z = B. 
Inserting (1551) into (BUI) , we find: 

g(B) = -— S m[np X (B)} (41) 

7T 1 — C 

Since the sign of g(B) is not allowed to change on the support, we must impose p = 1. 
The final answer for the resolvent can be written: 

W{zm = - i- (1 ~ c C)3/2 {- 7 =L= + cosh[7r( 7 /2 + ^/tf')]}> (42) 

and the corresponding pole density is: 



9(B) = ^-\ — - sm[7r X (B)} . (43) 

We can give an alternative definition for the function x(B) in terms of Jacobi elliptic functions [GW89 : 
if we set k' = vL — k 2 , 

x(B) = ±nd^(-^-). (44) 

The ratio k = i3 m i n /i? m ax is fixed as a function of c = tanh 2 (7T7/2) by (|36p with p = 1, and then we 
find from JSSJ) 

S max = — (45) 

7T7 1 — C 

To conclude, formula (|28|l delivers the expression of the front slope: 



^(-) = - 7 ^-n[| 7 sc-, 1 (-£-)], (46) 



where sc^., 1 is the reciprocal function of the Jacobi elliptic function scj.' = sn^ / cn^/ . These results are 
in agreement with the numerical resolution of , see Fig. [TJ 

3.3 Solution for -1 < c < 

As we explained, the only possibility in this case is -B m in = 0. It corresponds to the limit k — > 0, 
and thus K — > tt/2 in the former construction, in which the Jacobi elliptic function degenerates to 
trigonometric functions. We prefer to use the change of variable: 

i?:i:nv (47) 



cosh ip ' 



which, although not exactly the limit of (|30|) . is closely related. Thanks to (|47|) . the z-complex plane 
with its two cuts [— -B m ax,0] U [0, B max ] is mapped onto the strip T> = {0 < Im-0 < 7r} (see Fig. [3]). 
More precisely, the segment [0, B max ] ± iO is mapped to the half-line M=p ; the segment [— -B max , 0] ± iO 
is mapped to the half-line iir + M T ; the point z — oo is mapped to = iir/2. 
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17T — 00 





^■-plane 



-oo 




— B, 



oo 



A. 



17T - 





+0O 





Figure 3: Mapping between the z-plane with two merged cuts [— £? ma x,0] U [0, B max ] and the strip 
T> in the ^-plane (values in black). For special points in the ip-p\ane, we indicate below (in red) the 
corresponding value of z. The blue dashed (resp. green dotted) arrow indicates the image in the 
-0-plane of a path in the z-plane with positive (resp. negative) imaginary part. 



oj(ip) = W(z(ip)) defines a holomorphic function on the strip T>. Since the analysis now parallels 
the case c > 0, we go quickly over it. Continuity of the resolvent W(z) across [— -B ma x, 0] implies 



w(V0 = w(2i7r — ip) 

when Imi/) = tt, and (|26|) translates into 

cj(V>) + oj{— ip) + nuj(ip + in) = 1/v 



(48) 



(49) 



when tp belongs to R. These relations can be used to analytically continue ui(ip) in the whole complex 
plane, so that (|48|) - (j49"l) are satisfied everywhere, or equivalently (|48|) and 

ui(tp) +uj(tp + 2i7r) + nu(ip + in) = 1/v. (50) 

The general solution of (I4"51) - (l4l?l) can be written: 



1 



+ ff(V0 W +(V>) + g{-i>) w+(2i7r - ip), 



(51) 



v{2 + n) 

where g(tp) is an arbitrary in periodic, analytic function, and w+ is any analytic function satisfying: 

uj+(^ + m) = e 47r ^ + (?/>)■ (52) 
where we recall the parametrization (fT0|) n = — (e 17r7 + e~ 17r7 ). A possible choice is: 



U}+(lp) 



e (7-i)^ 
sinh( - 0) 



(53) 



Our change of variable (l4"Tl) sends the point z = at ip = ±oo, more precisely z ~ {B max /2) e T ^. 
So, we have to be careful that W(V0 does not grow too much when ip —> ±oo so that the pole density 
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g(B) remains integrable at B = 0. This condition is equivalent to demanding that W(ip) be o(e^) 
when ip — > +00. Being a in periodic function, g(tp) admits a Fourier expansion g(iji) = Spgz 5p e2p ^ ■ 
The integrability condition results in g p — whenever p ^ 0, 1. Besides, since the right-hand side of 
([4T)| is a constant, we find that uj(ip) is an entire function. We deduce that, in the solution we are 
looking for, g(ip) must be holomorphic, and such that the apparent pole of (|53l) at ip = m is absent 
in (|51[) . This gives a relation between go and g\ , and leads eventually to our solution in the form: 

= + + G cosh [ 7( ^ ~ ^ (54) 
for some constant G. The normalization condition (1271) becomes: 



iTV 

"WO = 0+- (ip - hr/2) + o(ip - m/2), (55) 

-t>max 

and provides two equations determining the two unknowns B max and G. We find that the maximal 
position of the poles is: 

4 , s 

Bmax = - t Nv. 56 

7 1 — c 

We thus find a unique solution for the resolvent: 

1 (1 - c) 3 / 2 / 1 \ 
W{Z m = ~v -c \JT=l ~ C ° Sh[7( ^ " i?r)] ) ' (57) 
and as expected, the corresponding pole density is positive: 

1 /I - <■ i . r . -1 (B 



0(B) = ± Ji— ^ - sinh Ucosir 1 (^)l , (58) 
2n V — c ^ L \ B /J 

The flame slope is retrieved from (|2"5|) : 

fjf> x (a:) = p= sinh 7sinh~ 1 f — _ (59) 



We observe a strong similarity with the results of the regime < c < 1, namely (|4"2"1) - (|4"6"|) . Going from 
c > to c < 0, one just need to replace c by — c, the elliptic functions by their trigonometric analog, 
and 7 is replaced by i7, so that trigonometric functions rather appears as hyperbolic functions. 

3.4 Discussion 

The first noteworthy point concerns — 1 < c < 0: the results (|56|) and ([58 ]) - ([59]) then coincide with 
those found in |JD12j . up to slightly different notations (e.g., 7 was called /1 whereby = s ' n ^ ) 

and thanks to such identities as 2sinh~ 1 ( sin 1 h g ) = sgn(£) In ( tanh ^/ 2 ) ) • 

If < c < 1, however, the above exact results do differ from those found in |JD12] by naive 
continuation from c < 0. In particular, (|4"o| and (l4Ul) - (l4~4l are free of the oscillations that crippled ([3]) 
and the corresponding density profile. To explain why the latter "solutions" were nevertheless so close 
to their numerical counterpart over the x + \B = O(Nv) range when < c < 0.6, one first notices 
that the maximum pole altitude B max in (|45l) and the quantity B* = — > vN featured in are 
related by 

% = — = l + 0(fc 2 ) (60) 

B* TT k^O 

As K' — ln(4/fc) + o(l) when k — > 0, 7 ~ must be small in (|36p for p = 1, which implies: 



4ex P (-^=) (61) 
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The ratio k = g mi " is thus exceedingly small if < c <C 1, and it remains small if < c < 0.6 thanks 
to the 7r 2 constant in (1611) : £>» and -B max are then nearly equal. Similarly, the difference between 
(|4*6"f and (O is strongly suppressed by 0(k 2 ) factors wherever \x\ ^> B mln ; this is enough for §5§ 
to accurately capture the crest slope profile, the exact maximum max|0 a; (x)| = cT 1 ! 2 (reached at 
\x\/B m i n ~ k^ 1 / 2 > 1) inclusive. 

The discarded unphysical profiles associated with integers p > 1 in ([36]) would have had an extra 
p in the exponent of (|61|) . And, as p — > oo at fixed < c < 1 in ([36jl leads to k — > and -B m in — > 0, 
the front slope <fi x (x) given by (|4"6")l would resume ©• Indeed, the only "error" in ([3]) was the implicit 
assumption that B min / B max — 0, as when c < 0. 

Next, for k 1", we have if' - vr/2 and if = ± In (tj^) + o(l), while tt7 = In ( j4j) + o(l) if 
c — » 1~ ; equations (j3"6")l and (|4"5l) then reduce to: 

(l - ~ 4^1^, B max ~ ^ (62) 

V -Bmax / 1-C 

quite in line with Appendix [A] Using the integral representation of nd^", 1 |WW 96, p. 494] in (|4"3")) - (l4*4")) 
and (f6"Tj) . one can show that (|4"3")) resumes the Marcenko-Pastur law when c — > 1~. 



4 Bicoalesced periodic fronts 

4.1 Integral equations for pole densities 

We now consider steady periodic cells, for which <fi(t, x) = —Vt+<f>(x) ; the drift velocity was calculated 
in jJou91j in terms of the neutral-to-actual wavelength v G]0, 1[ and the total number TV of pairs of 
poles involved V = 2Nv( ^-~Nv) ^ Tjp 0n re scaling of the space coordinate, we may assume that (j) x {x) 
is 27r-periodic. We will address more specifically the issue of bicoalesced periodic fronts. These are 
obtained as periodic solutions of the ZT equation which admit a pole decomposition, with two vertical 
piles of poles separated by a half-wavelength: 

Ni v_ N 2 v_ 

<j) x {x)= Y /" C R ; + Y , 1 ' C ., : (63) 

A large number of these bicoalesced solutions have been obtained in |Den06] in the Michelson- 
Sivashinsky case c = 0, all of which have been found to be stable with Neumann boundary conditions 
in the same paper. A solution of this type, obtained by direct numerical simulation, is presented in 
TBLOfl. Periodic configurations featuring a single pile of pole can be extracted from (|63|) by assum- 
ing N2 = 0, or assuming Ni = N2 and bk = Bk, and doubling the wavelength. The analysis of such 
solutions is parallel to that of isolated crests. Plugging (|63|) in the ZT equation fixes the prefactor 
equal to and yields coupled equations for the position of the poles. In the first pile, we have for 
any k = 1, .. . , Ni, 

Nl 1 n Nl 1 

gtanh(^)"2g tann (M) (64) 



Etanh(^^)-^gtanh(^±^) = /, 



N 2 „ N- 

n 

2 

m—1 rn—1 
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and symmetrically in the second pile, we have for any m = 1, . . . , iV 2 : 

N 2 N 2 



gtanh(^) tgtanh(^) (65) 

N i , D iVi l , d 

■Et-h(^^)-^tanh(^±^) = /, 

fc=l k=l 



where the constant in the right-hand side reads: 



-fl+ -^—Nu), N^N 1 +N 2 . (66) 
v V 1 — c / 



We shall study ([M|) - (p5)) for large crests, i.e. in the limit Ni,N 2 — > 00. We expect that the 
poles B\, . . . , Bn x of the first pile (resp. 61, ... , &at 2 of the second pile) get condensed on a segment 
■h = [B min , -Bmax] (resp. J 2 = [b m i n ,b max ]) witn a continuous nonnegative density Qi{B) (resp. £> 2 (&)), 
as discussed in Section [2] The densities are normalized as: 

/ Qi(B)dB = N u f g 2 (b)db = N 2 . (67) 

J J! JJ 2 

and the flame slope is retrieved as: 

, , , -2v r y giCgQdg' y g2 (fr')dfr' i 

^ = ~ R ^ A MSI) + i j2 tan(^)J (68) 
The relations (|64l) - ([65l) then turn into coupled integral equations for the densities: for all B £ Ji, 

Ql (B>)dB> n f gi (B')dB' 

- f 

J j 2 cotanh ( ) 2 J j 2 cotanh ( ^ ) 

and for all b € J2: 

g 2 (6')dfe' n y 2 (h')db' 



Jit anh(^) 2y, 7l tanh(™ 
02 (6')d& / n y g a (6')d6 / 



Ja tanh ( ^ ) 2 7 j2 tanh ( ^ ) 
Ql {B')dB' n y ^(B') dB' 



Jt cotanh(^) 2 7^ C otanh(^) 



(70) 
= /• 



Even though N — > 00, the full expression of / given in (|66|) is to be retained in (J69J)- (JTOJ) . The 
reason why steady periodic fronts need an at least 0(N) actual-to- neutral wavelength ratio 1/v will 
be discussed in § 14.61 



4.2 Reformulation by resolvents 

To take advantage of the periodicity, we perform a first change of variable: 

m tanh(B/2) tanh(fe/2) 
tanh( J B max /2) ' tanh(6 max /2) ' 

and define the resolvents for a complex variable z: 



(71) 



W)= ,«, W J*»t. (72) 

Z — 1 I r_ Z T 
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where J\ = [T m j n , 1] and J 2 — [r m i n , 1] are the image of the supports after (fTTj). The normalization 
conditions (|67p are now rephrased as: 

N = f 2tanh(E max /2) gl (5(T))dr = f dz 2 tanh(E max /2) W^z) 
1 J A 1 - tanh 2 (£ max /2)T 2 fj 1 2iir 1 - tanh 2 (B max /2)z 2 ' 

and the contour of integration can be moved at infinity to pick up residues at z = ±1. Therefore: 

N ^ = W i{ tanh(B max /2) ) ~ Wl ( tanh(B max /2) ) ' (?4) 

The same computation holds for N 2 , with Wi and £? max replaced in (I7il) by W2 and 6 ma x- Starting 
from the essential properties of resolvents emphasized in (|!M)) - (|2"5)) , the relations ([6"!))) - (fT0"| can be 
transformed after a tedious algebra into a Riemann-Hilbcrt problem for Wf. for any T 6 J±, 

W X {T + iO) +Wi(T — iO) + nWi(-T) + 2W 2 (1/PT) + nW 2 (-l/PT) = /, (75) 
and a symmetric equation for any r € J2, 

VK 2 (t + iO) + W 2 (t - iO) + nTy 2 (-r) + 2W 1 (1/Pt) +nff 2 (-l/Pr) = /, (76) 
where we have introduced the - yet unknown - constants: 

P = tanh(B max /2)tanh(6 max /2), (77) 

/ = /-t^-M, (78) 
1 — c 

M = f tanh 2 (B max /2)r gl (B(r))dT f tanh 2 (6 max /2) r g 2 (b(r))dr 

7,7, l-tanh 2 (B max /2)r 2 7j 2 1 - tanh 2 (6 max /2) r 2 ' ''' 

4.3 Symmetric case 

As of writing, we could not solve d75l) - d76l in their full generality. Yet, we did determine the solution 
when there is some symmetry between the two piles of poles, the same as for c = |JD08j . viz.: 

Ji = J 2 = [r, 1] and Qi(B(T)) = g 2 (b(T)) = g(T). (80) 

This means that Q\{B) and £> 2 (&) can be deduced from one another by the rescaling encoded in 
(fTI]) . We stress that in general Qi{B) ^ 92(B), unless £? m ax = &max, i.e. unless they have the 
same support (in which case (|68j) describes two copies of a monocoalesced periodic cell of length ir). 
Likewise, the number of poles in each pile need not be equal. When (pi)) is satisfied, we obviously have 
Wi{z) = Wi(z). The trick is now to introduce a new resolvent, which will incorporate simultaneously 
the terms involving z and 1/z in (|75[) - (I76I) . For this purpose, we perform a second change of variable: 

« = Y^Pf2> ( 81 ) 

already introduced in |JD08] to deal with the case c = 0, and which is invariant under T ■<-» 1/PT. 
Though not invertible on the whole complex T-plane, (|8ip sends bijectively the support [r, 1] of the 
density g in the T-plane, to the segment: 

r 1 

in the 77-plane. We can thus define: 

p{rj) = g(T( v )), (83) 
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which is a density supported on J, and next the resolvent for a complex variable z: 

W( Z ) = [ M*L. (84) 



A computation shows that it is related to ([721) by: 

W(y^) = W(z) + W(1/P Z ) + £ 2 7 + g jgf - (85) 

We deduce from (1751 (or (|76p) the relation, for all 77 G J: 

W(r) + iO) + W(r] - iO) + n W(-tj) = f. (86) 

The constant in the right hand-side combines (|77|) and the last term in and becomes quite simple 
at the end: 

f = f--^M, (87) 

1 — c 



where 

?2 



«2^_^(*(i /S ) + *(-i /S )), 



M = 

S = tanh(5 max /2) + tanh(6 max /2), 



and we have evaluated M thanks to Cauchy residue formula. We now recognize an equation of the 
same type as (|26|) . hence solvable. The main difference with Section [2] is that the data of the number 
of poles Ni and N2 are encoded in a different way in the resolvent. We can rewrite ([74]) in terms of p 
only: 

f (S+ , )r^L, i = l,2, (89) 



N t = 

where e\ = +1, e 2 = — 1, and: 

D = tanh( J B max /2) - tanh(& max /2) (90) 
So, the total number of poles is retrieved in a simple way from W: 

N = W(1/S)- W(-l/S), (91) 
whereas the difference between the two piles is only given implicitly: 



Nx-N 2 f D p{rj)d.n 



(92) 



y/1 - 4Pr, 2 (1- SW)' 

Eventually, the flame slope (|68l) can be expressed directly in terms of the resolvent: if we first define 



s(x ) = tan^/2) 

tanh(£ max /2)-tanh(& max /2)tan 2 (x/2)' 

we find after some algebra that: 

*M--£/,i^-i^M*<W.>>]- W 

Before coming to the determination of <p x (x), let us recapitulate the unknowns and the parameters for 
those symmetric solutions. Initially, the supports [B m - m , B max ] and [6 m i„, 6 max ] of the pole densities 
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in the first and second pile were unknown. The symmetric case we are studying amounts to consider 
supports with the same aspect ratio, namely: 

^ _ tanh(5 min /2) _ tanh(b min /2) 
tanh(S max /2) tanh(6 max /2) ' 

We are left with two unknowns, which can be encoded in the variables P and S defined in (|7T[) and 
([88]) . On the other hand, c (or equivalcntly n or 7), and the number of poles in each pile Ni and N2, 
count as free parameters. 

4.4 Solution for < c < 1 

If we assume -B m in, b m in > 0, we have r] m i a > (see (|82"T) ). We can repeat the analysis of Section [321 
We first perform the change of variable 

, , Vmin (l+P)r 

Z = Vmin sn k ^, k = = 2 , (96) 

'/max x ~r ± I 

and w(V') = W(z(ip)) can be extended as a holomorphic function of ip in the whole complex plane. 
We deduce that ?7 m i n > implies c > 0, and the quantization condition: 

jK' — 2pK for some positive integer p. (97) 

We already know that the positivity of the pole density requires p = 1. Besides, the general solution 
of ([55} in terms of 10 is obtained from where 1/v (which appeared in the right-hand side of &o§) 
is replaced by / given in (|66)) and (|78p. Exploiting the parametrization ([9]) of n in terms of c, we 
arrive at: 

u{il>) = ^ + i[w(^i / s)+w(-^i / s)]+GcoBh[7r(7/2 + V'/«' / )] ) (98) 

where G is a constant to determine and tpi/s = i-^' + ^ with < S < K is the unique point satisfying 
(IM)) for z = 1/S 1 . This can be rephrased as: 

snfe ^ = ITp v 2 J' ( ^ 

We have a priori five unknowns: S*, P, r as discussed in Section |4~3[ the constant G, and ^ [uj(tp 1 / s ) + 
uj(— 4>i/s)] (which is in fact irrelevant in the final result for the density and the flame slope) appearing 
in (198[) . We already have two normalization conditions (liTT1) - (l9"2l and the quantization condition ([9"7| . 
A fourth equation comes from the fact that W(rj = 00) = (see ([84)) ). hence uj(iK') — 0, and fixes 
the value of ^[ui(ipi/g) + w(— ipi/s)] in terms of the other unknowns. The fifth and last equation is 
obtained by consistency of (|M)l : if we specialize to ip = ±ipi/s and consider the half-sum of the two 
equations, we find: 

-Gcosh(7r 7 /2)cosh(^) = -^L-. (100) 

Besides, (EHI) can be rewritten: 

- 2Gsinh(7T7/2) sinh (^j = N. (101) 

Solving for G and 4>i/s m the system (|100|) - (I101[) . and taking into account the expression ([6"6"| for /, 
yields: 

F (1 - c) 3/2 I 4c 

G =-Tu c ■ F = Jl+—Nv(l-Nv), (102) 
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— theoretical solution 
■- ■- numerical solution 




x coordinate 



Figure 4: Bicoalesced periodic front shape <f>{x) with c > 0. The plot assumes (Ni,N 2 ) = (200, 100), 
1/v = 600.5 and c = 0.5, We display the comparison between the numerical resolution of ([M]) - ([rj5|) 
and integration of the theoretical formula (|106[) . 



and: 



1 + P 



sn fc 



K' 



tanh 1 



2NvJc 



1 - c(l - 2A^) 



(103) 



Then, the system of equations (l92l) . ([97|) . (|103[) determines implicitly S, P and r in terms of N\,Ni and 
c, and it can be solved at least numerically. The final answer for the reduced density of poles (|80|) , (f83| 
reads: 



1 1 — c F 



2tt 



sin[7rx(?7)L 



where we have introduced: 

x{v) = jp nd^, 1 ( 
And from (|94|) . we deduce the flame slope: 

F 



1 + Pr 2 



7) , ife' = Vl-k 2 



4>x{x) 



Sill 



^(l + Pr 2 



s(x) 



(104) 



(105) 



(106) 



where s(x) and F were defined in (|93p and (|102[) . A flame shape <f>(x) corresponding to (| 106[) is plotted 
in Fig. |H 



4.5 Solution for -1 < c < 

In this regime, we must have -B m i n = 6 m i n = 0, i.e. rj m m = or equivalently r = 0, and we can repeat 
the analysis of Section 13.31 We first perform the change of variable: 

?7max 



cosh lp ' 



(107) 



the function uj(ip) = W(z(^)) can be extended as a holomorphic function in the whole 'i/'-plane, and 
has the general expression (|54[) where 1/v is replaced - as in Section [4.4| - by the appropriate constant 
term: 

w ^ = irrz + \n^/s)+u(;™ -i>x/s)] +Gcosh[7(^-i7r)], (108) 
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J i I i I I I I I i L 

-3-2-10 12 
x coordinate 



Figure 5: Bicoalesced periodic front shape <f>{x) with c < 0. The plot assumes (Ni,N 2 ) — (200, 100), 
\jv = 600.5 and c = —0.25. We display the comparison between the numerical resolution of (|64 p -(|65 p 
and integration of the theoretical formula (|113|) . 



where G is some constant to determine, and ipi/s = K 71 /^ ~ with < i5 < 7r/2 is the unique point 
satisfying (| 107[) for z = 1/S. In other words: 



sin fo) = = tanh 

v ' 1 + P 



( ^max "T "max \ 

V 2 )' 



(109) 



Since the discussion to determine the unknowns is also very similar to Section 14. 4| we only give the 
results: 

F (1-c) 3 / 2 



G= , 

4z/ -c 

with the same constant F appearing in (I102p . and: 



5' 



rl 



1 + P 



■ tan 



2Niv 



I - r(] - 2-Vl/)) 



(110) 



(111) 



Together with the implicit relation (1921) . Eq. (jllip fixes S and P, and thus the solution of our problem. 
The final result for the reduced density of poles (|80 ]) . ([83P reads: 



1 1-c F 



2ir V — c v 
and the corresponding flame slope is: 

F 



sinh 



7 cosh 



1 



(1 + P)r? 



<t>x(x) 



sinh 



7 sinh 1 



1 



(1 + P)«(x) 



(112) 



(113) 



where s(x) and F were defined in (|77|) and p02l) . As expected, the resulting expression for 4> x (x) is 
symmetric under (a;, P ma x) (x — tt, i) max ). A flame shape (j>(x) corresponding to f| 1 13[) is plotted in 
Fig. [S] 



4.6 Discussion 

As is patent on comparing ([TT2 ]) -([TT3 ]) with pIM)) - ([TUS)) . using -1 < c < instead of < c < 1 
again merely amounts to selecting the relevant expression of 7, and to replacing the elliptic functions 
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I 1 

theoretical solution c=0.25 
numerical solution c=0.25 
theoretical solution c=0.5 
numerical solution c=0.5 




x coordinate 



Figure 6: Mono-coalesced periodic front shapes <f>{x) with c > 0. The plots assume (Ni,N 2 ) — (100, 0) 
and X/v = 199.5. For the two values c = 0.25 (top) and c = 0.5 (bottom), we show the comparison 
between the numerical resolution of (j64|) - (|65|) and integration of the theoretical formula (| 1 06[) . 



by suitable hyperbolic ones, in the pole density and flame slope. Both situations thus share several 
trends, which are appropriate to discuss first. 

Since siik(K) = 1 = sin(7r/2), £? max + 6 max — > oo follows in the limit 2vN — > 1, from (|M|) and 
(fT03)) if < c < 1, or (fT09|) and (|TTT]) if -1 < c < 0. Next, (J92J) indicates that only the uppermost end 
of the most populated pile (namely, -B max if N\ > ^2) actually goes to infinity in either case ; since 
B nra cannot go further vertically, steady bicoalesced periodic patterns cannot support an arbitrary 
large number of poles per wavelength [Jou91j . just like when c = [TFH85, JD08]. This can be traced 
back to (I64p - (I55|) . for example with &at 2 < Bn 1 — > 00 if N% > N2, in which limit the condition 



N = N x + N 2 = N opt (u) = L- 



1 



2u 



(114) 



is obtained. Note that N Qpt (2v) < N opt (i')/2, but the difference is asymptotically small for large 
wrinkles, as 2vN ov ^{y) ~ 1 for v —> + . Conversely, a large actual-to- neutral wavelength ratio 
l/v> (2N — 1) is required for a total population of 2N bicoalesced poles to stay steady "in" a cell. 



The analysis of Section @] of the large steady wrinkles governed by the ZT equation (fT]) contain as 
limiting cases all the results available so far. Monocoalesced 27r-periodic cells correspond to 2z/A^2 — > 
at fixed vNi (see Fig [6]). In the present continuous approximation of the pole equations, & ma x then 
vanishes and this amounts to replacing: 



1 



(1 + P)r) 
1 

(l + P)s(x) 



tanh(5 max /2) 
tanh(B/2) ' 
tanh(5 max /2) 
tan(cc/2) 



in (I112p and (|113|) respectively. For such two-crested periodic patterns, tan ^jff x m , 

B 



/2) 



(115) 



plays the part 



did for isolated crests, see P5 ]) -(|1S |1 and O-O. Also, if -1 < c < 0, ([TlT1)- pT?|) then resume 
the expressions found in [JD08J . The isolated crests constitute a further degeneration of (I115P and 
formally corresponds to 2vN\ ^ 1 and restricting oneself to a region of the complex plane where 
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LB I = 0{N x v) < 1. Then £ max = 0{N x v) < 1 by (gg) or p9]) . 

1 B roax tanh(-B roax /2) S n 



(1 + P)r? ' S ' tan(a;/2) ' a: ' ^ 116 ^ 

in (1115[) . the prefactor F defined in (|102l) goes to 1, and the density and the flame slope reduce to 



(g3D-(|46D or d58D-<[59j. 

For equally populated piles of poles (i.e. Ni = N 2 ), we have b max = i? max , whence (|81"T) and (|93l) 
simplify to: 

1 _ tanh 5 max 1 _ tanh B maK 

(l + P)s(x)~ tanx ' (1 + P)t]~ tanhS ' *- ' 

The comparison with (I115P shows that (| 1 1 7[) describes two copies of a monocoalesced cell of wavelength 
7r: this was actually expected from the outset in view of the discrete pole equations (l64|) - (l65l) with 
Ni = N 2 , since tanh(a) + - — r = t — ™— r . 

1 z ' V / 1 tanh(a) tanh(2a) 

The main qualitative difference between the c < and c > cases undoubtedly is about the 
tips of large front wrinkles. Even when N\$ — > oo (and keeping vN\, 2 = 0(1)), these are smooth 
maxima of 4>(x) if c > 0, with 0(-B m in) radii of curvature that quickly increase with c if < c -C 1 
(see Fig. EE]), and are sharp cusps otherwise (see Fig. [5]). Specifically, the centered tips locally have 
(j) x {x) oc — sgn(x) \x\~ a and thus 

4>{x) - (f)(0) oc -\x\ 1 ~ a , a = -tanh" 1 ( % /^). (118) 

with < a < 1/2 for —1 < c < 0. The MS value c = is marginal and was already known 
to give logarithmic cusps [TFH85| IJD08 . While affecting the crests (— H[<f> x ] > 0), the value of c 
also modifies the front troughs (minima of 4>(x), where H[4> x ] > 0), but differently. This can be 
qualitatively understood from (fTJ), by transferring the nonlocal nonlinearity cH[<fi x ] 2 to the right-hand 
side : if c < 0, the resulting "effective driving term" of instability, — H[(j) x ](l + %U[<f>x]), is reduced 
(or enhanced) at troughs (or crests) compared to its MS counterpart alone and tends to render them 
rounder (or sharper) ; c > has opposite influences. 

And, importantly again, whatever the value of — 1 < c < 1 is, the analytical predictions of Section|4] 
perfectly agree with the numerical resolutions of (|64l) - (|rj5|) for the discrete poles Bf, and b m , and with 
the bicoalesced front slopes ensuing from (|63[) and numerical quadrature. In particular, <f>(x) is found 
for c > to have 4 inflexion points whenever Nx/N 2 ^ 0, oo (see Fig. @|, and only 2 inflexion points 
otherwise ; in either case, max |0 x (a;)| = ^= is the same at all of them and only depends on c > and 
Nv. 



5 Conclusion 

We have investigated with complex analysis methods solutions of the ZT equation with free parameter 
— 1 < c < 1, describing stationary flame shapes in the limit of large wrinkles. We rederived the known 
expression (|59| for a front with an isolated crest when — 1 < c < 0, and obtained its counterpart (|46|) 
in the regime < c < 1. We also obtained new expressions describing bicoalesced, space-periodic 
fronts with some symmetry (see § I4.3[) in the regime < c < 1 (|106[) or — 1 < c < (|113[) , which 
are in agreement with numerical simulations (resp. Fig. [4] and Fig. [5]). The results obtained in the 
limit case of monocoalesced periodic fronts are also new (see Section 14.61 and Fig [6] for comparison to 
numerics). 

The preceding analyses do not exhaust all the theoretical problems as to (fTJ), for it admits even 
more general solutions than the symmetric bicoalesced ones studied in Section SJ For example, the 
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so-called interpolating solutions [GM031 IDen06] are still awaiting for detailed descriptions. Their 
existence can be inferred on noticing that, besides 2N bicoalesced poles vertically aligned at x = 
or i = 7T, iVoo extra pairs may stay in equilibrium at Zk — ±ioo + f or k — 1, . . . , if 

DO 

< 2i/N = 1 — vNoc < 1 ; with 2vN being further reduced, the Noo remote poles will move and be 
located at finite Im Zk |GM031 IDen06j . Although such solutions are likely unstable (as in the c = 
case), their computation in the limit v — > + and still uNx,2,oo = 0(1) with similar methods is under 
investigation. As another example, space-periodic configurations with more than 2 unequal piles per 
cell certainly exist if v is small: a single pair of poles may already stay in equilibrium near the trough 
of a base front if wide enough, thereby creating an incipient extra crest there [JD08] . Whether this is 
within reach of resolvent approaches when the extra poles get many constitutes a challenge still to be 
met, as some of the poles would again condense on curved arcs in such configurations. 

In the continuous description of, say, isolated crests, v could have been absorbed in the normal- 
ization of the pole density, since the flame slopes (|59p or (141)1) only depend on v and N through the 
combination Af = vN . This parameter Af could then be any positive real number, labeling a contin- 
uum of solutions. Without invoking pole decompositions, can one retrieve the selection of a discrete 
set of solutions (e.g. requiring N = M jv to be an integer) if small curvature effects (0 < v <C 1) 
are restored afterward ? This was answered in the affirmative for viscous fingering or needle-crystals 
[PelOOj : quantization then resulted from a solvability condition on short-scale steady front-shape per- 
turbations (local wavenumbers of order The corresponding procedure for flames has so far not 
been provided, and would be tantamount to performing directly a WKB analysis in (|20j) . 

A related matter concerns the solutions with integer p > 1 first encountered in S I3.2l but discarded 
to ensure positivity of the pole density. We have shown in § 12. 31 that they nonetheless lead to solutions 
of the "inviscid ZT equation", i.e. ([T]) without the v4> X x term. The p > 1 oscillating <fi x (x) profiles 
obtained for < c < 1 could then be viewed as steady nonlinear perturbations of the p = 1 solution 
(|46| , in a way reminiscent of what happens in needle-crystal growth [PelOOj . 

Besides, other equations than (JT|) may reduce to the inviscid ZT form in the large wrinkle limit, 
at least outside the cusps they might have. In case the same outer profiles as described by (|58|) - ([59| 
would apply, one would need to determine the appropriate Af that enters in -B max in such equations 
as (|56j) . One may guess that the inner front tip structure yields Af (e.g., by matching), and one may 
wonder whether its value is quantized. 

We also claim (the details are found in Appendix [Cj that the ZT equation in the form (|2"0"| appears 
in the /3-deformation of the O(n) random matrix model (|128l) . in the limit of large matrices. One may 
wonder whether this very matrix model for finite N (and its generalization where the eigenvalue lives 
in a certain region of the complex plane) could provide itself a statistical-mechanical description of a 
combustion problem. 



A The c — > 1 limit 

The then strong coupling between Bk of unlike signs in ([5]) makes the pole population split in two 
groups, separated by a distance of about 2h > v i±| that will exceed the width of each group for 
c-)T. With Bj K Q m —h to leading order, (J6)) written for B\, ... , Bjy simplifies to: 

A 2v „ 1 + c uN 



^ (B k +h)- (Bj +h) l-cB k + h 
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According to Stieltjes (see |Sze39] . the review [MMFMG07] and references therein), the solutions 
£,k = Bk + h to (|119p are the zeros of the associated Laguerre polynomial C^\C) , with a = — 1 + 2N . 
One may select h so that b-(-b ) ~ "5T7T ^ e exac t aX B = h, whereby £ = SL ^- = N y±| ; 

whatever N is, Bm — B\ oc J^_ c |IL92] is asymptotically smaller than h ~ when c — > 1~, as 
anticipated. If next N — > oo, the zeros of £^(£) are distributed according to a Marcenko-Pastur 
density |MP67j 



p(T,) CX V(Cmax ^min)^ (12Q) 



with £ max / m i n = iV(2 + a ± 2\/l + a) and a = linijv^oo 4& « 2 The ensuing pole density is 
nonzero only for i? m i n < < B ma:x , with -B max ( m ; n ) ~ ^Cmax(min) — ft > 0. 

B Origin of the quantization ( |36| ) 

We first recall the definition of the first Jacobi theta function: 

This series is absolutely convergent whenever the parameter r has positive imaginary part, it defines 
an entire function of w € C with the following properties: 

tfiO + l|r) = -0i(V>|t), (122) 
^(tfl + rlr) = -e-^^+^iHr), (123) 
= Q|r) = 0. (124) 

Therefore, we may build a function u; + satisfying (|35|) with a ratio of theta functions. The choice of 
uj + is arbitrary: indeed, the ratio of two functions satisfying (|35p is 2K- and 2ii4T'-periodic, so that it 
can be absorbed in the choice of g(ip) in (f34|) . In this appendix, we take: 

= ' T = ^- (125) 

This function has a simple pole at ip = and a simple zero at ?p = —^K' (modulo the translations by 
2K and 2iiT), unless: 

jK' = 2Kp for some integer p. (126) 



We are now in position to explain why the case considered in § 13.21 imposes the quantization 
condition (|126l) . Since the right-hand side in ([2^)1 was a constant, we were looking for a entire 
function to(ip). The latter was decomposed as (l34l) . hence: 

+^ ' e 2T7-l e 2 ^-12 + n V ; 

is also an entire function (we recall the assumption c ^ 0, thus 7 > 0). When (jl26[) is not satisfied, the 
choice of (|125j) would imply that 5(1/') has a simple pole at ip = —jK', and no other singularity in C. 
But no such function exists, because the total order of poles (modulo translations) of a meromorphic 
doubly-periodic function is at least 2. The quantization (|126[) was thus necessary for our problem to 
admit solutions. 
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C 0(n) matrix models and ZT equation 



Let us consider the statistical ensemble of N particles at position B\, . . . , Bn on the positive real axis, 
distributed according to the measure: 

J" Y[^ N iB l + B 3 y^ ( 128 ) 

where V(B) is a given, smooth function. For (3 — 2, this is the 0(n) model introduced by [GK89, Kos89 
in relation with the problem of counting configurations of self-avoiding loops on random discrete 
surfaces. This model has been studied intensively since then |KS92[ IEK 95, EK96, Borllj . and also 
appears in the context of quantum entanglement BN12 . When (3 = 2 and n is an integer, (|128[) is 
the measure induced on eigenvalues B±, . . . , Bn of a random hermitian, positive definite matrix B, 
and coupled to n other hermitian matrices Ai, . . . , A n with a joint distribution: 

dB f[ dA j e -^[v(B)+E? =1 ba»] _ (129) 

The quantity A(z) 2 = Yl i<: j \z% — zj\ 2 is characteristic of the eigenvalue distribution of random hermi- 
tian matrices whose full distribution is invariant under conjugation by a unitary matrix: it expresses 
the repulsion between eigenvalues of such a matrix taken at random. The product in the denominator 
of (|128[) arises from the integration over the matrices Aj (which have a Gaussian distribution accord- 
ing to (|129p ). In random matrix theory, (jl 29[) is called "the O(n) matrix model", and replacing A(z) 2 
by A(z)P is the "/3-deformation" . 

In this model, the B/~ are random. We can define again a resolvent as an expectation value against 
the measure (|128p : 

N , 

W(z) = (^—) t (130) 

i— 1 

We claim that, when N is large and with the choice V(B) — Bjv, W(z) satisfies an equation like (|2"0)l , 
which amounts to saying that 4> x (x) = yz~ Im[W(ix)] is solution of a ZT equation (see § 12. 3[) . 
It is convenient to introduce the two-points resolvent: 

W{Z 1 ,Z 2 ) = (Y J 7T 5-)- (131) 

\ zi - Bj z 2 - B k l 

In general W(zi,z 2 ) ^ W(z 1 )W(z 2 ). However, when N becomes large, the random B k for k = 
X, . . . ,N are distributed in a deterministic way, with some density g(B) supported on C. Hence: 

W(zi,z 2 ) ~ W{zx)W{z 2 ) wheniV->oo. (132) 

Using integration by parts, one can derive the Schwinger- Dyson relations valid for any N [Borllj : 

W(z, z) + W{-z, -z) + nW(z, -z) + (| - l) [W'(z) + W'(-z)] - - \W{z) + W{-z)] = 0. (133) 
This can be simplified when N is large owing to (|132[) : 

W(z) 2 + W(-z) 2 + nW(z)W(-z) + (| - l) [W'(z) + W'(-z)] - ~ [W(z) + W(-z)] = 0, (134) 
which coincides, up to the rescaling z —> 2 /p-i > w ^ n d2Q]). 
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